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Numerical simulations of the approach to the singularity in vacuum spacetimes are presented 
here. The spacetimes examined have no symmetries and can be regarded as representing the general 
behavior of singularities. It is found that the singularity is spacelike and that as it is approached, 
the spacetime dynamics becomes local and oscillatory. 
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A longstanding problem in general relativity has been 
to find the general behavior of singularities. Several re- 
sults, both analytical 1] and numerical 2] have been ob- 
tained. However, for most of the results, the space- 
times have one or more symmetries. Only for scalar 
field matter have results been found when there are no 
symmetries. Q For vacuum, Belinski, Lifschitz and 
Khalatnikov 5] (BKL) have conjectured that the generic 
singularity is local, spacelike and oscillatory. This con- 
jecture has been reformulated and put more precisely by 
Uggla et a/.0 

In this paper are presented numerical simulations of 
the approach to the singularity in vacuum spacetimes 
with no symmetries. The results support the BKL con- 
jecture. 

The system evolved here is essentially that of 
reference[6J but specialized to the vacuum case and with 
a slightly different choice of gauge. Here the spacetime 
is described in terms of a coordinate system (t, x 1 ) and a 
tetrad (eo,e a ) where both the spatial coordinate index 
i and the spatial tetrad index a go from 1 to 3. Choose 
eo to be hypersurface orthogonal with the relation be- 
tween tetrad and coordinates of the form eo — N~ 1 dt 
and e a = e Q ! 9, where N is the lapse and the shift is 
chosen to be zero. Choose the spatial frame {e Q } to be 
Fermi propagated along the integral curves of e . The 
commutators of the tetrad components are decomposed 
as follows: 

[e ,e Q ] = ii a e - (HSa 13 + a a l3 )e p (1) 
[e a ,ep] = (2a[ Q (5 (3 ] 7 + e a psn 5l )e 1 (2) 

where n a/3 is symmetric, and cr al3 is symmetric and trace 
free. 

Scale invariant variables are defined as follows: 
{d Q ,d a } = {e ,e a }/H 

{E a \Z aP , A a , N a0 ] ee {e a \ a af3 , a a ,n af3 }/H (3) 

q + 1 ee — do In H and r a = —d a In H. 

Finally choose the lapse to be N = H -1 . The relation 
between scale invariant frame derivatives and coordinate 



derivatives is 8q — d t and d a — E a l di. From the vacuum 
Einstein equations one obtains the following evolution 
equations: 



d t E a l = FjE p l 

d t r a = F a p rp + d a q 

d t A a = F a fj A f3 + ±d /5 £ Q/3 

<9 f S Q/3 = (q-2)Z afj -2N <a 1 N /3> ' 1 + N^N <a 



(4) 
(5) 
(6) 



/3> 



+ d <a rP> - 8 <a A 0> + 2r <a A f}> 

+ e^ s<a {d 1 -2A 1 )N p> s (7) 

d t N al3 = qN afi + 2^ a s N^ s - e^d^s (8) 

d t q = [2(q - 2) + | (2A a - r a )8 a - \d a d a ] q 

- ld a r a + lA a r a + lr p d a Z al3 

- 2Z al3 W af3 (9) 

Here angle brackets denote the symmetric trace- free part, 
and F a p and W a p are given by 



(10) 



a/3 



— 2 



N ai Nf - \N\N aP + \d a A 



ld a rp-\ei 5 a {d 1 -2A 1 )N fj5 (11) 
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In addition to the evolution equations, the variables sat- 
isfy constraint equations as follows: 

= (C com )^ ee 2(d [a - r [a - A [a )E p] l 

- e a p 5 N^Ej (12) 
= C G ee 1 + \(2d a - 2r a - 3A a )A a - \N a pN af} 

+ ^(iV Q a) 2 -|^S a3 (13) 

= (C c ) a = dpY 1 a(} + 2r a -EV' 3 ~ 3^S Q/3 

- e a ^N ps ^ 5 (14) 
= C q ee q - i£^E a/3 + \d a r a - \A a r a (15) 
= (Cj)° ee {dp - rp)(N a ? + e a ^A,) 

- 2ApN af3 (16) 
= (C W ) Q = [e a ^{d p - A p ) - A^]r 7 (17) 

We want a class of initial data satisfying these con- 
straints that is general enough for our purposes but 
simple enough to find numerically. We use the York 
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methodQ and choose H to be constant, r a and N a p 
to vanish and the following form for the other vari- 
ables: Ej = (t/j- 2 /H)Sj,A a = 2c l d a <if> and E Q/3 = 
(— ^~ 6 /-ff)diag(£i, S2, E3). Then the constraints are 
satisfied provided 9 l S l = and q — ^ip~ 12 H~ 2 t, l t,i and 



d'diip = l(6H 2 ip 5 -E l E^~ 7 ) (18) 

We use the following solution for £^ 

Ei = a 2 cos y + (13 cos z + b 2 + b 3 
E 2 = o>i cos x — a 3 cos z + b 1 — 63 
E 3 = — a! cosx — a 2 cosy — — 6 2 (19) 

where the a; and 6^ are constants. We consider space- 
times with topology T 3 x R with each spatial slice hav- 
ing topology T 3 . In terms of the coordinates we have 
< x < 2n with and 2ir identified (and correspond- 
ingly for y and z). 

The numerical method used is as follows: each spatial 
direction corresponds to n + 2 grid points with spacing 
dx = 2ir/n. The variables on grid points 2 to n + 1 are 
evolved using the evolution equations, while at points 1 
and n+2 periodic boundary conditions are imposed. The 
initial data is determined once equation l|18|l is solved. 
This is done using the conjugate gradient method. The 
evolution proceeds using equations il with the excep- 
tion that the term (5 — 2q)C q is added to the right hand 
side of equation © to prevent the growth of constraint 
violating modes. Spatial derivatives are evaluated us- 
ing centered differences, and the evolution is done us- 
ing a three step iterated Crank-Nicholson method^] (a 
type of predictor-corrector method). In equation @ the 
highest spatial derivative term is —^d a d a q which gives 
this equation the form of a diffusion equation. Note that 
diffusion equations can only be evolved in one direction 
in time, in this case the negative direction which corre- 
sponds to the approach to the singularity. Stability of 
numerical evolution of diffusion equations generally re- 
quires a time step proportional to the sqare of the spatial 
step. However, the constant of proportionality depends 
on the coefficient of the second spatial derivative. To en- 
sure stability, we define -E max to be the maximum value 
of \E a l \ (over all space and over all a and 1) and then 
define dti = — j(dx / E ma ^) and dt 2 = — ■gdx. The time 
step dt is then chosen to be whichever of dt\ and dt 2 has 
the smaller magnitude. 

Before presenting numerical results, it is helpful to con- 
sider what behavior to expect as the singularity is ap- 
proached (that is as t — > —00). First denote the eigen- 
values of E Q ^ by (Ei, E 2 , E3). Then suppose that at suf- 
ficiently early times the time averages of q — E^ are all 
positive. Then the time averages of the eigenvalues of 
F a p are all positive. Since we are evolving in the neg- 
ative time direction, this should lead (through equation 
(QJ) to an exponential decrease in E a l . However, since 
all spatial derivatives appear in the equations through 
d a = E a l di we would expect the spatial derivatives to 



become negligible. That is, at each spatial point the 
dynamics becomes that of a spatially homogeneous cos- 
mology: the approach to the singularity is local. Note 
that this does not mean that the spacetime is becoming 
homogeneous. Spatial variation is not becoming small; 
however since all spatial derivatives appear in the evolu- 
tion equations through E a l di and since E a l is becoming 
small, the effect of the spatial derivatives on the evolution 
is becoming negligible, the positivity of the time aver- 
ages of the eigenvalues of F a p should also lead (through 
equations ijoKil) ) to exponential decrease in r a and A a . 
Thus we expect that as the singularity is approached, the 
dynamics is described by a much simpler version of evo- 
lution equations 1 111)1) and constraint equations (|I2I17|) 
where r a and A a and all spatial derivatives are dropped. 

This simpler set of equations describes homogeneous 
cosmologies, and has been treated extensively in the liter- 
ature on such spacetimes. [H ITo| Here we simply summa- 
rize the important attributes. First, from equation (|I4|I it 
follows that Ti a p and N a p commute and therefore have 
a common basis of eigenvectors. From equations <|7I8[) 
it follows that this basis (which we will call the asymp- 
totic frame) is not changed under time evolution. One 
can therefore express equations (|7I8|I in the asymptotic 
frame yielding equations for the E^ (eigenvalues of £ Q /j) 
and the Ni (eigenvalues of N a p). During a time period 
where all the N are negligible, it follows that all the 
Ej are constant. Such a time period is called a Kasner 
epoch. During a Kasner epoch, two of the Ni are decay- 
ing and one is growing. The growing eigenvalue of N a p 
eventually gives rise to a transition (called a "bounce" ) to 
another Kasner epoch. It is helpful to define a quantity 
u by 



E^E^-yE^ 



G 



8ki 2 (l 



(1 



u 2 Y 



(20) 



(there is a unique u > 1 provided the quantity on the 
left hand side of the equation is between -6 and 6). The 
quantity u is constant in each Kasner epoch and changes 
from epoch to epoch as follows: u ^ u — 1 if u > 2 and 
u — > l/(u— 1) if 1 < u < 2. This rule is called the u 
map. 

We now turn to the numerical simulations. All runs 
were done in double precision on a SunBlade 2000 with 
11 = 50 (except for an examination of resolution which 
also used n — 25). The equations were evolved from t = 
to t = —130. The initial value of H was -| corresponding 
to an initial trace of extrinsic curvature equal to — 1 . The 
hi were chosen with 63 = and b\ and b 2 given so that 
the spacetime would be a Kasner spacetime with u = 2.3 
if the a,i vanished. For the runs presented here the 
were given as (0.2, 0.1, 0.04). 

We would like to know whether E^, r a and A a become 
negligible as the singularity is approached. In figurenare 
plotted the maximum values (over all space, a and i) of 
In I-Eq, 1 !, In \r a \ and In \ A a \ as functions of time. Note the 
steep decrease in these quantities near t ~ —20. This 
indicates that r a , A a and the spatial derivatives become 
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FIG. 1: maximum values of In \E a l \ (solid line), In \r a \ (dot- 
ted line) and In \A a \ (dot-dashed line) vs time 



FIG. 3: components of E a/ 3 vs time, in the asymptotic frame: 
Ei (solid line), E2 (dotted line) and E3 (dot-dashed line) 
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FIG. 2: Constraint C q vs time, for n — 50 (solid line) and 
n — 25 (dot-dashed line) 



FIG. 4: components of N a p vs time, in the asymptotic frame: 
iVi (solid line), N2 (dotted line) and N3 (dot-dashed line) 



negligible for t < —20. (the failure of the quantities plot- 
ted in figure 1 to continue to decrease is most likely due 
to unresolved small scale spatial structure to be discussed 
below) . Thus the interesting part of the dynamics can be 
seen by looking at the development of the variables at a 
single point as a function of time. We now present data 
of that form. The behavior at the spatial point chosen is 
typical. 

The behavior of a constraint is presented in figure |21 
Here what is plotted is In \C q \ at the spatial point as a 
function of time. The solid line is for the n = 50 run 
and the dot-dashed line is for the n = 25 run. The finer 
resolution yields a smaller value for the constraint; but 
the resolution is not good enough to be in the convergent 
regime. Note that the time dependence of the two runs 
becomes increasingly out of sync. This is due to the 
chaotic nature of homogenous cosmologies. Also note 
that the shape of the constraint for the two runs does not 
completely match. This may be due to unresolved small 
scale structure to be discussed below. Similar results 
were obtained for the other constraints. 

Figures |3 and 0] show respectively the diagonal compo- 
nents of E Q( g and N a p in the asymptotic frame. Though 



not shown here, I have also examined the off-diagonal 
components of both E a( g and N a p in this frame. For 
times t < — 20 these off-diagonal components become 
negligible. This demonstrates that for these times the 
asymptotic frame is essentially unchanged and the quan- 
tities in figures [3] and 0] are eigenvalues of S"^ and N a p 
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FIG. 5: u vs time 
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respectively. Note that for t < — 20 the behavior of the 
components of £ a( g consists of epochs where they are con- 
stant punctuated by short bounces where they change 
rapidly. Note too that for t <J — 20 the behavior of the 
components of N a p is that they are negligible during the 
epochs of constant £ a /3 and that one component of N a f) 
becomes non- negligible at each bounce. This is exactly 
what we would expect from the approximation of using 
the equations for homogeneous spacetimes. 

We now turn to the behavior of the quantity u. In 
figure [S] is plotted u as a function of time. Note that u 
undergoes a series of bounces when the components of 
£ a do. The sequence of values of u begining at t ~ —20 
is 1.279,3.583,2.584,1.584,1.712. This sequence obeys 
the u map. 

In summary, these simulations provide strong support 
for the BKL conjecture. The initial data evolved have no 
symmetry and can be regarded as generic. The evolution 
shows that spatial derivatives become negligible and the 
time dependence goes over to the well studied behavior 
of a general homogeneous cosmology. This dynamics is 
oscillatory consisting of a series of Kasner epochs punc- 
tuated by short bounces. The details of the oscillations 
given by the behavior of the quantity u are in agreement 
with what would be expected for locally homogeneous 
spacetimes. 

We now consider what remains to be done on this sub- 
ject. First recall that bounces occur when one of the JVj 
grows exponentially. If that iVj initially vanishes on a 
surface S then we would expect bounces on either side 
of S, but not on S itself. This gives rise to a small scale 
structure which can be seen (crudely) in the results of 
these simulations, and which has been well studied in the 
Gowdy spacetimes. [Til Il2^ A corresponding study for the 
case of no symmetry will require high resolution and is 
work in progress. 



The present work only treats vacuum spacetimes. 
There is a general expectation that for most types of 
matter (a scalar field is an exception) the influence of the 
matter on the dynamics should become negligible as the 
singularity is approached. The formalism of reference 
is for a fluid with equation of state P = kp for constant k. 
Thus a fairly straightforward generalization of the simu- 
lations reported here would be to do simulations of the 
approach to the singularity for P = kp perfect fluid. 

Another question is whether there are any residual ef- 
fects of the spatial derivatives as the singularity is ap- 
proached. Comparison of the full evolution to one where 
the spatial derivatives are set to zero after a time to yields 
differences: the sequence of bounces is the same, but they 
occur more rapidly in the full evolution. This may be due 
to the spatial derivatives' increasing the values of the iV, 
and thus hastening the time when the JVj become large 
enough to cause a bounce. 

Finally, note that this simulation is for a spatially 
closed universe. Since the result is that the dynamics 
become local as the singularity is approached, these re- 
sults should describe at least a portion of the singularity 
in any generic spacetime, including asymptotically flat 
spacetimes that undergo gravitational collapse to form 
black holes. Nonetheless, there is a body of work[H Q 
that indicates that when a black hole forms, that portion 
of the singularity that is near the event horizon is null 
(or asymptotically null). It would be good to extend the 
methods presented here so that they are able to treat 
collapse in asymptotically flat spacetimes and the near 
horizon properties of the singularities formed. 

I would like to thank Mark Miller, Beverly Berger, 
Woei-Chet Lim, John Wainwright, Lars Andersson, Jim 
Isenberg and G. Comer Duncan for helpful discussions. 
This work was partially supported by NSF grant PHY- 
0244683 to Oakland University. 



[1] for a review see A. Rendall, "Theorems on existence and 
global dynamics for the Einstein equations" Living Re- 
views in Relativity (2002-6) 

[2] for a review see B. Berger, "Numerical Approaches to 
spacetime singularities" Living Reviews in Relativity 
(2002-1) 

[3] L. Andersson and A. Rendall, Commun. Math. Phys. 

218, 479 (2001) 
[4] D. Garfinkle, Phys. Rev. D65, 044029 (2002) 
[5] V. Belinskii, I. Khalatnikov and E. Lifschitz, Sov. Phys. 

Usp. 13, 745 (1971) 
[6] C. Uggla, H. van Elst, J. Wainwright and G.F.R. Ellis, 

Phys. Rev. D68, 103502 (2003) 
[7] J.W. York, Phys. Rev. Lett. 26, 1656 (1971) 
[8] W. Press, S. Teukolsky, W. Vetterling and B. Flannery, 



Numerical Recipes in FORTRAN, second edition (Cam- 
bridge University Press, Cambridge, 1992) 
[9] M. Choptuik, in Deterministic Chaos in General Relativ- 
ity, edited by D. Hobill, A. Burd and A. Coley (Plenum, 
New York, 1994), pp. 155-175 
[10] "Dynamical systems in cosmology" Edited by J. Wain- 
wright and G.F.R. Ellis, Cambridge University Press, 
1997 

[11] B. Berger and D. Garfinkle, Phys. Rev. D57, 4767 (1998) 
[12] A. Rendall and M. Weaver, Class. Quantum Grav. 18, 
2959 (2001) 

[13] E. Poisson and W. Israel, Phys. Rev. D41, 1796 (1990) 
[14] A. Ori and E. Flanagan, Phys. Rev. D53, 1754 (1996) 



